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ABSTRACT 

We use techniques from nonparametric function estimation theory to extract the density pro- 
files, and their derivatives, from a set of iV-body dark matter halos. We consider halos generated 
from ACDM simulations of gravitational clustering, as well as isolated, spherical collapses. The 
logarithmic density slopes 7 = d\ogp/d\ogr of the ACDM halos are found to vary as power-laws 
in radius, reaching values of 7 « —1 at the innermost resolved radii, ~ 10 r„i r . This behavior is 
significantly different from that of broken power-law models like the NFW profile, but similar to 
that of models like de Vaucouleurs'. Accordingly, we compare the Af-body density profiles with 
various parametric models to find which provide the best fit. We consider an NFW- like model 
with arbitrary inner slope; Dehnen & McLaughlin's anisotropic model; Einasto's model (identi- 
cal in functional form to Sersic's model but fit to the space density); and the density model of 
Prugniel & Simien that was designed to match the deprojected form of Sersic's R}l n law. Over- 
all, the best-fitting model to the ACDM halos is Einasto's, although the Prugnicl-Simicn and 
Dehnen-McLaughlin models also perform well. With regard to the spherical collapse halos, both 
the Prugniel-Simien and Einasto models describe the density profiles well, with an rms scatter 
some four times smaller than that obtained with either the NFW-like model or the 3-parameter 
Dehnen-McLaughlin model. Finally, we confirm recent claims of a systematic variation in profile 
shape with halo mass. 

Subject headings: dark matter — galaxies: halos — methods: N-body simulations 
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1. Introduction 

A fundamental question is the distribution of 
matter in bound systems (galaxies, galaxy clus- 
ters, dark matter halos) that form in an expand- 
ing universe. Early work on the self-similar col- 
lapse of (spherical) primordial overdensitics re- 
sulted in virializcd structures having density pro- 
files described by a single power law (e.g., Fill- 
more & Goldreich 1984; Bertschinger 1985; Hoff- 
man 1988). Some of the first TV-body simulations 
were simple cold collapse calculations like these 
(e.g. van Albada 1961; Aarseth 1963; Henon 1964; 
Peebles 1970). It was quickly realized that, given 
appropriately low but non-zero levels of initial ran- 
dom velocity, the end state of such systems de- 
parted from a simple power law, resembling in- 
stead the de Vaucoulcurs (1948) R 1 ^ 4 profiles ob- 
served in elliptical galaxies (e.g. van Albada 1982; 
Aguilar & Merritt 1990). As iV-body techniques 
improved, the logarithmic profile slopes of cold 
dark matter (CDM) halos, simulated in hierarchi- 
cal merger models, were also observed to steepen 
with increasing radius (e.g., West, Dekel & Oem- 
ler 1987; Frenk et al. 1988; Efstathiou et al. 1988). 
Dubinski & Carlberg (1991) adopted Hernquist's 
(1990) double power law model (itself a modifica- 
tion of Jaffe's 1983 model) to describe these den- 
sity profiles. This empirical model has an inner 
logarithmic slope of —1 and an outer logarithmic 
slope of —4. It was introduced as an analytical 
approximation to the deprojected form of de Vau- 
couleurs' (1948) profile. Navarro, Frenk, & White 
(1995) modified this to give the so-called NFW 
model that has an outer logarithmic slope of —3 
rather than -4, while Moore et al. (1998, 1999) 
suggested that a further variation having an inner 
logarithmic slope of —1.4 or —1.5 may be more 
appropriate. 

The density profiles of TV-body halos typically 
span only ~ 2 decades in radius, between the virial 
radius and an inner limit set by the TV-body reso- 
lution. It has long been clear that other functional 
forms might fit such limited data as well or bet- 
ter than the NFW or Moore profiles. Recently, 
Navarro et al. (2004) argued for a model, like de 
Vaucouleurs', in which the logarithmic slope varies 
continuously with radius: 



dlnp 
d\nr 




i.e. 

p(r) oc exp(-ylr a ) , (2) 

where r_2 is the radius at which the logarithmic 
slope of the density is —2 and a is a parameter 
describing the degree of curvature of the profile. 
Merritt et al. (2005) pointed out that this is the 
same relation between slope and radius that de- 
fines Sersic's (1963, 1968) law, with the differ- 
ence that Sersic's law is traditionally applied to 
the projected (surface) densities of galaxies, not to 
the space density. Merritt et al. further showed 
that the (space) density profiles of a sample of TV- 
body halos were equally well fit by equation (2), 
or by a deprojected Sersic profile, and that both 
of these models provide better fits than an NFW- 
like, double power-law model with a variable inner 
slope. Hence, Sersic's law - the function that is so 
successful at describing the luminosity profiles of 
early-type galaxies and bulges (e.g. Caon, Capac- 
cioli, & D'Onofrio 1993; Graham & Guzman 2003, 
and references therein), and the projected density 
of hot gas in galaxy clusters (Demarco et al. 2003) 
- is also an excellent description of TV-body halos. 

(To limit confusion, we will henceforth refer to 
equation (2) as "Einasto's r 1 /™ model" when ap- 
plied to space density profiles, and as "Sersic's 
i? 1 /™ model" when applied to projected density 
profiles, with R the radius on the plane of the sky. 
The former name acknowledges Einasto's (1965, 
1968, 1969) early and extensive use of equation (2) 
to model the light and mass distributions of galax- 
ies (see also Einasto & Haud 1989). In addition, 
we henceforth replace the exponent a by 1/n in 
keeping with the usage established by Sersic and 
de Vaucouleurs.) 

In this paper, we continue the analysis of alter- 
natives to the NFW and Moore profiles, using a 
new set of TV-body halos. Among the various mod- 
els that we consider is the Prugniel-Simien (1997) 
law, first developed as an analytic approximation 
to the deprojected form of the Sersic R}l n profile. 
Apart from the work of Lima Neto et al. (1999), 
Pignatelli & Galletta (1999), and Marquez et al. 
(2000, 2001), the Prugniel-Simien model has re- 
ceived little attention to date. Demarco et al. 
(2003) have however applied it to the gas den- 
sity profiles of 24 galaxy clusters observed with 
ROSAT, and Terzic & Graham (2005) showed that 
it provides a superior description of the density 
profiles of real elliptical galaxies compared with 
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either the Jaffe or Hernquist models. As far as 
we are aware, ours is the first application of the 
Prugniel-Simien model to A-body halos. 

As in Merritt et al. (2005), we base our model 
evaluations on nonparametric representations of 
the A-body density profiles. Such representations 
are "optimum" in terms of their bias-variance 
tradeoff, but are also notable for their flexibil- 
ity. Not only do they constitute (i) "stand-alone," 
smooth and continuous representations of the den- 
sity and its slope: they are also well suited to (ii) 
inferring best-fit values for the fitting parameters 
of parametric functions, and (iii) comparing the 
goodness-of-fit of different parametric models, via 
the relative values of the integrated square error or 
a similar statistic. The more standard technique 
of computing binned densities is suitable (though 
inferior) for (ii) and (iii) but not for (i), since the 
density is given only at a discrete set of points 
and the derivatives poorly defined; while tech- 
niques like Sarazin's (1980) maximum-likelihood 
algorithm provide a (perhaps) more direct route 
to (ii) but are not appropriate for (i) or (iii). Re- 
cently, nonparametric function estimation meth- 
ods have been applied to many other problems in 
astrophysics, including reconstruction of the CMB 
fluctuation spectrum (Miller et al. 2002), dynam- 
ics of dwarf galaxies (Wang et al. 2005), and re- 
construction of dark matter distributions via grav- 
itational lensing (Abdclsalam, Saha & Williams 
1998). Application of nonparametric methods to 
the halo density profile problem is perhaps over- 
due, especially given the importance of determin- 
ing the inner density slope (Diemand et al. 2005). 

In §2 we introduce the data sets to be analyzed. 
These consist of A-body simulations of ten ACDM 
halos and two halos formed by monolithic (nearly 
spherical) collapse. (Moore ct al. 1999 have dis- 
cussed the similarity between the end state of cold 
collapse simulations and hierarchical CDM mod- 
els.) In §3 we present the nonparametric method 
used to construct the density profiles and their 
logarithmic slopes. §4 presents four, 3-parameter 
models, and §5 reports how well these empirical 
models perform. Our findings are summarized in 
Section 7. 

In Paper II of this series (Graham et al. 2006a), 
we explore the Einasto and Prugniel-Simien mod- 
els in more detail. Specifically, we explore the 
logarithmic slope of these models and compare 



the results with observations of real galaxies. We 
also present the models' circular velocity pro- 
files and their p/cr 3 profiles. Helpful expressions 
for the concentration and assorted scale radii: 

the radius where 
the circular velocity profile has its maximum value 
— are also derived. Because the Prugniel-Simien 
model yields the same parameters as those com- 
ing from Scrsic-model fits, we are able to show in 
Paper III (Graham et al. 2006b) the location of 
our dark matter halos on the Kormcndy diagram 
(/x e vs. log -R e ), along with real galaxies. We addi- 
tionally show in Paper III the location of our dark 
matter halos and real galaxies and clusters in a 
new log(pe) — l°g(i? c ) diagram. 

2. Data: Dark matter halos 

We use a sample of relaxed, dark matter halos 
from Diemand, Moore, & Stadel (2004a,b). De- 
tails about the simulations, convergence tests, and 
an estimate of the converged scales can be found 
in those papers. Briefly, the sample consists of 
six, cluster-sized halos (models: A09, B09, C09, 
D12, E09, and F09) resolved with 5 to 25 mil- 
lion particles within the virial radius, and four, 
galaxy-sized halos (models: GOO, G01, G02, and 
G03) resolved with 2 to 4 million particles. The 
innermost resolved radii are 0.3% to 0.8% of the 
virial radius, r v ; r . The outermost data point is 
roughly at the virial radius, which is defined in 
such a way that the mean density within r v ; r 
is 178rj^ 45 /9 cri t = 98.4p crit (Eke, Cole, & Frenk 
1996) using Vt m = 0.268 (Spergel et al. 2003). The 
virial radius thus encloses an overdensity which is 
368 times denser than the mean matter density. 
We adopted the same estimates of the halo centers 
as in the Diemand ct al. papers; these were com- 
puted using SKID (Stadel 2001), a kernel-based 
routine. 

In an effort to study the similarities between 
cold, collisionless collapse halos and CDM halos, 
we performed two additional simulations. We dis- 
tributed 10 7 particles with an initial density pro- 
file p{r) oc r 1 , within a unit radius sphere with 
total mass 1 (Mil) and 0.1 (M35). The parti- 
cles have zero kinetic energy and the gravitational 
softening was set to 0.001. Each system collapsed 
and underwent a radial-orbit instability (Merritt 
& Aguilar 1985) which resulted in a virialized, 
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triaxial/prolate structure. The lower mass halo, 
M35, collapsed less violently over a longer period 
of time. 

3. Nonparametric estimation of density 
profiles and their derivatives 

Density profiles of TV-body halos are commonly 
constructed by counting particles in bins. While 
a binned histogram is a bona-fide, non-parametric 
estimate of the "true" density profile, it has many 
undesirable properties, e.g. it is discontinuous, 
and it depends sensitively on the chosen size and 
location of the bins (see, e.g., Stepanas & Saha 
1995). A better approach is to view the particle 
positions as a random sample drawn from some 
unknown, smooth density p(r), and to use tech- 
niques from nonparametric function estimation to 
construct an estimate p of p (e.g., Scott 1992). In 
the limit that the "sample size" N tends to infin- 
ity, such an estimate exactly reproduce the den- 
sity function from which the data were drawn, as 
well as many properties of that function, e.g. its 
derivatives (Silverman 1986). 

We used a kernel-based algorithm for estimat- 
ing p(r), similar to the algorithms described in 
Merritt & Tremblay (1994) and Merritt (1996). 
The starting point is an estimate of the 3D den- 
sity obtained by replacing each particle at position 
Vi by a kernel of width hi, and summing the kernel 
densities: 



P(r) 



N 

2^ h* 



i=i 



(3) 



Here m, is the mass associated with the ith par- 
ticle and K is a normalized kernel function, i.e. a 
density function with unit volume. We adopted 
the Gaussian kernel, 
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The density estimate of equation (3) has no im- 
posed symmetries. We now suppose that p(r) = 
p(r), i.e. that the underlying density is spherically 
symmetric about the origin. In order that the 
density estimate have this property, we assume 
that each particle is smeared uniformly around 
the surface of the sphere whose radius is rj. The 



spherically-symmetrized density estimate is 

p(r) = £™±JvjMtoOK($>) 

(5b) 
(5c) 



|r - Ti\ 

rf + r 2 — 2rri cos 9 



where 9 is defined (arbitrarily) from the r^-axis. 
This may be expressed in terms of the angle- 
averaged kernel K, 

i r /2 

K(r,n,hi) = — d<f> 

4?r J-ir/2 

x d6 sin0 K (hi^^Jrf + r 2 - 2rr~cos^a) 



rf + r A 



2rnp , (6b) 



as 



N 



i=l a i 

Substituting for the Gaussian kernel, we find 

-l 



K(r,n,hi) = 



1 



(2tt)3/2 \hl 
x cxp [-{r 2 + r 2 )/2h 2 } sinh(rr,/^). ( 8 ) 

A computationally preferable form is 

-l 



K = 



1 / rn 



x {exp [-(n - r) 2 /2h 2 ] - exp [-(r t + r) 2 /2h 2 ] } (9) 

Equations (7) and (9) define the density estimate. 
Typically, one sets up a grid in radius and eval- 
uates p(r) discretely on the grid. However we 
stress that the density estimate itself is a contin- 
uous function and is defined independently of any 
grid. 

Given a sample of N positions and particle 
masses drawn randomly from some (unknown) 
p(r), the goal is to construct an estimate p(r) that 
is as close as possible, in some sense, to p{r). In the 
scheme just described, one has the freedom to ad- 
just the N kernel widths hi in order to achieve this. 
In general, if the hi are too small, the density es- 
timate will be "noisy," i.e. p(r) will exhibit a large 
variance with respect to the true density; while 
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if the hi are too large, the density estimate will 
be over-smoothed, i.e. there will be a large bias. 
(Of course the same is true for binned histograms, 
although in general the bias-variance tradeoff for 
histograms is less good than for kernel estimates.) 
If the true p(r) were known a priori, one could ad- 
just the hi so as to minimize (say) the mean square 
deviation between p(r) and p(r). Since p(r) is not 
known a priori for our halos, some algorithm must 
be adopted for choosing the hi. We followed the 
standard practice (e.g. Silverman 1986, p. 101) of 
varying the hi as a power of the local density: 

hi = h [ppiiot(n) / ffP" , (10) 

where p p i\ t(r) is a "pilot" estimate of p(r), and 
g is the geometric mean of the pilot densities at 
the Since the pilot estimate is used only for 
assigning the hi, it need not be diffcrcntiable, and 
we constructed it using a nearest-neighbor scheme. 

The final density estimate p(r) is then a func- 
tion of two quantities: ho and a. Figure 1 illus- 
trates the dependence of p(r) on h when the ker- 
nel algorithm is applied to a random sample of 10 6 
equal-mass particles generated from an Einasto 
density profile with n = 5, corresponding to typi- 
cal values observed in Merritt et al. (2005). Each 
of the density profile estimates of Figure 1 used 
a = 0.3. As expected, for small ho, the estimate 
of p(r) is noisy, but faithful in an average way to 
the true profile; while for large ho, p(r) is a smooth 
function but is biased at small radii due to the av- 
eraging effect of the kernel. For a = 0.3, The 
"optimum" ho for this sample is ~ 0.05 r c , where 
r c is the half-mass radius coming from the Einasto 
model (see Section 4.2). 

In what follows, we will compare the nonpar a- 
metric estimates p(r) derived from the iV-body 
models with various parametric fitting functions, 
in order to find the best-fitting parameters of the 
latter by minimizing the rms residuals between the 
two profiles. For this purpose, any of the den- 
sity estimates in Figure 1 would yield similar re- 
sults, excepting perhaps the density estimate in 
the uppermost panel which is clearly biased at 
small radii. In addition, we will also wish to char- 
acterize the rms value of the deviation between 
the "true" profile and the best-fitting parametric 
models. Here it is useful for the kernel widths to 
be chosen such that the residuals are dominated by 
the systematic differences between the parametric 



o 



0.1 






I ' 

0.03 






1 1 1 

0.01 


< 




0.003 " 




' I 1 1 

0.001 



0.01 0.1 1 



'A. 

Fig. 1. — Nonparametric, bias- variance trade- 
off in the estimation of p(r) using a single sam- 
ple of 10 6 radii generated from a halo having 
an Einasto r 1 /™ density profile with n = 5 
(see Section 4.2). From top to bottom, h = 
(0.1, 0.03, 0.01, 0.003, 0.001) r c ; all estimates used 
a = 0.3 (see equations 7, 9, and 10). 
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and nonparametric profiles, and not by noise in 
p(r) resulting from overly-small kernels. We veri- 
fied that this condition was easily satisfied for all 
of the A-body models analyzed here: there was 
always found to be a wide range of (ho, a) values 
such that the residuals between p(r) and the para- 
metric function were nearly constant with vary- 
ing radial coordinate. This is a consequence of 
the large particle numbers (> 10 6 ) in the A-body 
models, which imply a low variance even for small 
ho- 

As discussed above, quantities like the deriva- 
tive of the density can also be computed directly 
from p(r). Figure 2 shows nonparametric esti- 
mates of the slope, d\ogp/d\ogr, for the same 
10 6 particle data set as in Figure 1. We computed 
derivatives simply by numerically differentiating 
p(r); alternatively, we could have differentiated 
equation (9). Figure 2 shows that as ho is in- 
creased, the variance in the estimated slope drops, 
and for ho « 0.2 r c the estimate is very close to the 
true function. We note that the optimal choice of 
ho when estimating derivatives is larger than when 
estimating p(r) (~ 0.2 r c vs. ~ 0.05 r c ); this is a 
well-known consequence of the increase in "noise" 
associated with differentiation. Figure 2 also illus- 
trates the important point that there is no need 
to impose an additional level of smoothing when 
computing the density derivatives (as was done, 
e.g., in Reed et al. 2005); it is sufficient to increase 
h - 

3.1. Application to the A-body halos 

Figure 3 shows, using a — 0.3 and ho — 0.05 r e 
(left panel) and a = 0.4 and ho = 0.05 r c (right 
panel) 1 , the nonparametric estimates of p(r) (left 
panel) and j(r) = d\ogp/d\ogr (right panel) for 
the ten A-body halos. Figure 4 shows the same 
quantities for the two data sets generated from 
cold collapses. We stress that these plots - espe- 
cially, the derivative plots - could not have been 
made from tables of binned particle numbers. For 
most profiles, the slope is a rather continuous func- 
tion of radius and does not appear to reach any 
obvious, asymptotic, central value by ~ 0.01r v ; r . 
Instead, j(r) varies approximately as a power of 
r, i.e. log 7 vs. logr is approximately a straight 



x Wc have purposely used a small value of ho to avoid any 
possibility of biasing the slope estimates. 
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Fig. 2. — Five estimates of the logarithmic slope 
of an Einasto r 1 /" halo, derived via differentia- 
tion of p(r). The same sample of 10 6 radii was 
used as in Figure 1. From top to bottom, ho = 
(0.3, 0.2, 0.1, 0.05, 0.03) r ; each estimate used a = 
0.4 (see equation 10). Dashed lines show the true 
slope. 
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line. Accordingly, we have fitted straight lines, via 
a least-squares minimization, to the logarithmic 
profile slopes in the right-hand panels of Figures 3 
and 4. The regression coefficients, i.e. slopes, are 
inset in each panel. (These slope estimates should 
be seen as indicative only; they are superseded 
by the model fits discussed below.) In passing 
we note that such a power-law dependence of 7 
on r is characteristic of the Einasto model, with 
the logarithmic slope equal to the exponent l/n. 
Noise and probable (small) deviations from a per- 
fect Einasto r 1 /" model are expected to produce 
slightly different exponents when we fit the den- 
sity profiles in the following Section with Einasto's 
r 1 /" model and a number of other empirical func- 
tions. 

The slope at the innermost resolved radius is al- 
ways close to —1, which is also the slope at r = in 
the NFW model. However there is no indication in 
Figure 3 that j(r) is flattening at small radii, i.e., 
it is natural to conclude that iV-body simulations 
of higher resolution would exhibit smaller inner 
slopes. On average, the slope at r v i r is around —3, 
but there are large fluctuations and some halos 
reach a value of —4, as previously noted in Die- 
mand et al. (2004b). The reason for these fluctu- 
ation may be because the outer parts are dynami- 
cally very young (i.e. measured in local dynamical 
times) and they have only partially completed the 
violent relaxation into a stable, stationary equilib- 
rium configuration. We are not able to say with 
any confidence what the slopes do beyond r v ; r . 

4. Empirical models 

In this section, we present four parametric 
density models, each having three independent 
parameters: two "scaling" parameters and one 
"shape" parameter. We measured the quality of 
each parametric model's fit to the nonparamet- 
ric /3(r)'s using a standard metric, the integrated 
square deviation, 



d(logr) [logp(r) - logp param (r)] (11) 



where p pai a,m is understood to depend on the var- 
ious fitting parameters as well as on r. Equa- 
tion (11) is identical in form to the Cramer- von 
Mises statistic (e.g. Cox & Hinkley 1974, eq. 6), 
an alternative to the Kolmogorov-Smirnov statis- 



tic for comparing two (cumulative) distribution 
functions. 

We chose to evaluate this integral by discrete 
summation on a grid spaced uniformly in logr; 
our measure of goodness-of-fit (which was also 
the quantity that was minimized in determing the 
best-fit parameters) was 



T m S 2 



83 = log 10 [p(r,-)/ft> 



.fo)] 



(12a) 
(12b) 



with m = 300. With such a large value of m 
the results obtained by minimizing (12a) and (11) 
are indistinguishable. We note that the quantity 
A 2 in equation (12a) is reminiscent of the stan- 
dard x 2 , but the resemblance is superficial. For 
instance, A 2 as defined here is independent of to in 
the large-m limit (and our choice of to = 300 puts 
us effectively in this limit). Furthermore there is 
no binning involved in the computation of A 2 ; the 
grid is simply a numerical device used in the com- 
putation of (11). 

4.1. Double power-law models 

Hernquist (1990, his equation 43) presented a 5- 
parameter generalization of Jaffe's (1983) double 
power-law model. Sometimes referred to as the 
(a, (3, 7) model, it can be written as 



1+ - 



P{r) = p s 2 

(13) 

where p s is the density at the scale radius, r s , 
which marks the center of the transition region 
between the inner and outer power-laws having 
slopes of —7 and — (3, respectively. The parame- 
ter a controls the sharpness of the transition (see 
Zhao et al. 1996; Kravtsov et al. 1998; and equa- 
tions (37) and (40b) in Dchncn & McLaughlin 
2005). Setting (a,/3,7)=(l, 3, 1) yields the NFW 
model, while (1.5, 3, 1.5) gives the model in Moore 
et al. (1999). Other combinations have been used, 
for example, (1,3,1.5) was applied in Jing & Suto 
(2000) and (1,2.5,1) was used by Rasia, Tormen, 
& Moscardini (2004). 

In fitting dark matter halos, Klypin et al. (2001, 
their figure 8) have noted a certain degree of de- 
generacy when all of the 5 paramaters are allowed 
to vary. Graham et al. (2003, their figures 3 and 
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Fig. 3. — Nonparametric estimates of the density p(r) (left panel) and the slope d\ogp/d\ogr (right panel) 
for the ten TV-body halos of Table 1. The virial radius r V i r is marked with an arrow. Dashed lines in the 
right hand panels are linear fits of log(— d\ogp/d\ogr) to logr; regression coefficients are also given. 
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Fig. 4. — Nonparametric estimates of p(r) (left panel) and dlogp/dlogr (right panel) for the two "collapse" 
models. Dashed lines in the right hand panels are linear fits of log(— dlogp/dlogr) to logr. 



4) have also observed the parameters of this em- 
pirical model to be highly unstable when applied 
to (light) profiles having a continuously changing 
logarithmic slope. Under such circumstances, the 
parameters can be a strong function of the fitted 
radial extent, rather than reflecting the intrinsic 
physical properties of the profile under study. This 
was found to be the case when applied to the dark 
matter halos under study here. We have therefore 
chosen to constrain two of the model parameters, 
holding a fixed at 1 and (3 fixed at 3. 

In recent years, as the resolution in A-body 
simulations has improved, Moore and collabora- 
tors have found that the innermost (resolved) log- 
arithmic slope of dark matter halos has a range 
of values which are typically shallower than —1.5: 
recently obtaining a mean value (± standard de- 
viation) equal to —1.26 ± 0.17 at 1% of the virial 
radius (Diemand, Moore, & Stadel 2004b). At the 
same time, Navarro et al. (2004) report that the 
NFW model underestimates the density over the 
inner regions of most of their halos, which have 
innermost resolved slopes ranging from —1.6 to 
—0.95 (their Figure 3). A model with an outer 
slope of —3 and an inner slope of —7 might there- 
fore be more appropriate. Such a model has been 
used before and can be written as 



p(r) = 



2 3 "T/ 



{r/r s )i{l + r/r s f-i' 



(14) 



The total mass of this model is infinite however. 

We have applied the above (1, 3, 7) model to 
our dark matter density profiles, the results of 



which are shown in Figure 5 for the A-body halos, 
and in the upper panel of Figure 6 for the cold col- 
lapse models. The rms scatter, A, is inset in each 
figure and additionally reported in Table 1. 

4.1.1. Two-parameter models 

Recognizing that galaxies appear to have flat 
inner density profiles (e.g., Flores & Primack 1994; 
Moore 1994), Burkcrt (1995) cleverly introduced 
a density model having an inner slope of zero and 
an outer profile that decayed as r~ 3 . His model is 
given by the expression 



p(r) = 



Por s 



(r + r s )(r 2 + r s 2 ) 



(15) 



where p is the central density and r s is a scale 
radius. Application of this model in Figure 7 re- 
veals that, with only 2 free parameters, it does not 
provide as good a fit to the simulated dark mat- 
ter halos as the (1, 3, 7) model presented above. 
The hump-shaped residual profiles in Figure 7 sig- 
nify the model's inability to match the curvature 
of our density profiles. (It is important to point 
out that Burkert's model was introduced to fit the 
observed rotation curves in low surface brightness 
galaxies, after the contribution from the baryonic 
component had been subtracted out, a task which 
it performs well.) 

As noted previously, the NFW (a, (3, 7)=(1, 3, 
1) model also has only two parameters: p s and r s . 
Because this model is still often used, we apply it 
to our halos in Figure 8. Comparison with Figure 5 
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Table 1 
Three-Parameter Models 



Halo 


r s 


log p s 


7 


A 


r e 


log p e 


lEin 


A 


Rc 


log p' 


raps 


A 


Id. 


kpc 


M Q pc" 3 




dcx 


kpc 


M Q pc" 3 




dcx 


kpc 


M© pc" 3 




dex 






(1,3,7) 








Einasto 


r .l/n 






Prugnicl-Simien 














Cluster-sized halos 












A09 


626.9 


-3.87 


1.174 


0.025 


5962. 


-6.29 


6.007 


0.015 


2329. 


-2.73 


3.015 


0.021 


B09 


1164. 


-4.75 


1.304 


0.037 


17380. 


-7.66 


7.394 


0.041 


4730. 


-3.34 


3.473 


0.038 


C09 


241.8 


-3.27 


0.896 


0.040 


1247. 


-4.95 


3.870 


0.030 


738.9 


-2.55 


2.192 


0.016 


D12 


356.1 


-3.82 


1.251 


0.026 


2663. 


-6.02 


5.939 


0.020 


1232. 


-2.52 


3.147 


0.019 


E09 


382.5 


-3.96 


1.265 


0.033 


2611. 


-6.06 


5.801 


0.032 


1231. 


-2.62 


3.096 


0.030 


F09 


233.9 


-3.51 


1.012 


0.030 


1235. 


-5.26 


4.280 


0.025 


697.3 


-2.63 


2.400 


0.017 












Galaxy-sized halos 












GOO 


27.96 


-3.16 


1.163 


0.020 


189.0 


-5.22 


5.284 


0.023 


114.4 


-2.02 


3.135 


0.028 


G01 


35.34 


-3.36 


1.275 


0.029 


252.6 


-5.51 


5.873 


0.028 


146.0 


-2.01 


3.425 


0.032 


G02 


53.82 


-3.59 


1.229 


0.034 


391.4 


-5.74 


5.725 


0.031 


214.9 


-2.34 


3.243 


0.036 


G03 


54.11 


-3.70 


1.593 


0.028 


405.6 


-5.98 


7.791 


0.023 


229.1 


-1.47 


4.551 


0.024 












Spherical collapse halos 












Mil 


0.0175 


2.66 


0.006 


0.223 


0.244 


0.27 


3.426 


0.043 


0.187 


2.57 


2.445 


0.051 


M35 


0.0180 


1.62 


0.030 


0.249 


0.240 


-0.70 


3.214 


0.059 


0.185 


1.47 


2.301 


0.061 



Note. — Col.(l): Object Id. Col.(2)-(5) (1, 3, 7) model (equation 13 and 14) scale radius r s , scale density p s , inner profile slope 
7, and rms scatter of the fit. Col.(6)-(9) Einasto r 1 /™ model half-mass radius r c , associated density p c , profile shape riEin, and rms 
scatter of the fit. Col.(10)-(13) Prugnicl-Simien model scale radius R c , scale density p' (the spatial density p c at r = R e is such 
that p e = p'e~ b ), profile shape npg, and rms scatter of the fit. Note: The radius and density units do not apply to Mil and M35. 
For each halo, of the three models shown here the model having the lowest residual scatter is high-lighted in bold. 
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Fig. 5. — Residual profiles from application of the 
3-parameter (1, 3, 7) model (equation 14) to our 
ten, TV-body density profiles. The virial radius 
is marked with an arrow, and the rms residual 
(equation 12a) is inset with the residual profiles. 



reveals that the NFW model never performs better 
than the (1, 3, 7) model; the residuals are ^50% 
larger and sometimes twice as large. Importantly, 
the large-scale curvature observed in many of the 
NFW resdiual profiles (Figure 8) reveals that this 
model does not describe the majority of the halos, 
and that the (1, 3, 7) model should be preferred 
over the NFW model. 

An alternative 2-parameter expression has re- 
cently been studied by Dehncn & McLaughlin 
(2005, their equation (20b); see also Austin et al. 
2005). It is a special case of a more general family 
of models — which we test next — when the ve- 
locity ellipsoid at the halo center is isotropic and 
p/u r is a (special) power law in radius, varying 
as 7-- 35 / 18 . This 2-parameter density model is an 
(a, P, 7) = (4/9, 31 /9, 7/9) model given by 



PV ' ; (r/r a )V9[i + ( r / ra )4/9]6' ^ 

and is applied in Figure 9. It clearly provides a 
much better match to the dark matter halo den- 
sity profiles in comparison with the previous 2- 
parameter model over the fitted radial range, but 
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Table 2 

Three- Parameter Models (cont.) 



Halo 




logp s 


7' 


A 


Id. 


kpc 


Mq pc~ 3 




dcx 




Anisotropic Dchncn-McLaughlin 


(Eq.17) 






Cluster-sized 


halos 




A09 


722.7 


-2.21 


0.694 


0.013 


B09 


1722. 


-3.30 


0.880 


0.040 


C09 


207.0 


-1.34 


0.241 


0.047 


D12 


322.8 


-1.95 


0.683 


0.022 


E09 


330.4 


-2.04 


0.669 


0.034 


F09 


193.6 


-1.56 


0.350 


0.036 






Galaxy-sized 


halos 




GOO 


20.89 


-1.11 


0.422 


0.017 


G01 


25.88 


-1.28 


0.568 


0.023 


G02 


43.05 


-1.60 


0.581 


0.027 


G03 


30.20 


-1.34 


0.849 


0.024 




Spherical collapse halos 




Mil 


0.025 


4.23 


0.00 


0.179 


M35 


0.025 


3.21 


0.00 


0.206 



Note.— Col.(l): Object Id. Col.(2)-(5) Dehnen- 
McLaughlin (their equation 46b) scale radius r s , scale 
density p s , inner profile slope 7', and rms scatter of 
the fit. Note: The radius and density units do not 
apply to Mil and M35. When the rms scatter is lower 
than the value obtained with the other 3-parameter 
models, it is given in bold. 
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Table 3 
Two-Parameter Models 



Halo 


r a 


log po 


A 


r s 




A 




log/O s 


A 


Id. 


kpc 


M Q pc- 3 


dex 


kpc 


M pc- 3 


dex 


kpc 


Mq pc" 3 


dex 






Burkert 






NFW 




Isotropic Dehnen-McLaughlin (equation 16) 












Cluster-sized halos 








A09 


114.0 


-1.65 


0.242 


419.8 


-3.50 


0.042 


933.7 


-2.43 


0.018 


B09 


145.2 


-2.23 


0.247 


527.2 


-4.03 


0.068 


1180.0 


-2.97 


0.042 


C09 


96.16 


-1.74 


0.181 


284.4 


-3.42 


0.042 


554.3 


-2.27 


0.091 


D12 


68.39 


-1.62 


0.230 


213.3 


-3.34 


0.051 


409.1 


-2.17 


0.026 


E09 


77.09 


-1.80 


0.215 


227.0 


-3.46 


0.053 


428.2 


-2.28 


0.037 


F09 


80.17 


-1.85 


0.181 


229.0 


-3.49 


0.030 


438.2 


-2.32 


0.066 












Galaxy-sized halos 








GOO 


10.12 


-1.56 


0.139 


22.23 


-2.94 


0.024 


34.43 


-1.59 


0.037 


G01 


10.28 


-1.54 


0.152 


23.12 


-2.95 


0.038 


36.53 


-1.61 


0.031 


G02 


14.06 


-1.66 


0.183 


36.39 


-3.22 


0.044 


63.06 


-1.96 


0.035 


G03 


09.35 


-1.32 


0.179 


19.54 


-2.68 


0.066 


26.98 


-1.23 


0.025 










Spherical collap 


se halos 






Mil 


0.0261 


3.01 


0.203 


0.0309 


2.23 


0.233 


0.0234 


4.31 


0.244 


M35 


0.0265 


1.98 


0.231 


0.0314 


1.20 


0.259 


0.0236 


3.29 


0.269 



Note. — Col.(l): Object Id. Col.(2)-(4) Burkert (1995) model scale radius r s , central density po, and rms scatter 
of the fit (using m — 2 in the denominator of equation refEqChi). Col.(5)-(7) NFW (1, 3, 1) model scale radius r s , 
scale density po, and rms scatter of the fit (using m — 2). Col.(8)-(10) Dehnen-McLaughlin (2005, their equation 20b) 
model scale radius r s , associated density p s , and rms scatter (using m — 2). This model has an inner and outer, 
negative logarithmic slope of 7/9 w 0.78 and 31/9 w 3.44, respectively. Note: The above radius and density units do 
not apply to Mil and M35. For each halo, the 2-parameter model with the lowest residual scatter is high- lighted in 
bold. 
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Fig. 6. — Residual profiles from the application of 
seven different parametric models (see Section 4) 
to our "cold collapse" density halos, Mil and M35. 
Einasto's model, which has the same functional 
form as Sersic's model, is labelled 'Sersic' in this 
Figure. In the lower panel, the solid curve cor- 
responds to the 2-parameter model from Dchnen 
& McLaughlin (2005), and the dashed curve cor- 
responds to their 3-parameter model. The rms 
residual (equation 12a) is inset in each figure. 
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Fig. 7. — Residual profiles from application of 
Burkert's 2-parameter model (equation 15) to our 
dark matter density profiles. 
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Fig. 8. — Residual profiles from application of the 
2-parameter NFW (1,3,1) model to our dark mat- 
ter density profiles. 
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the rms scatter reveals that it does not perform as 
well as the (1, 3, 7) model, nor can it describe the 
'spherical collapse' halos (Figure 6). We therefore, 
in the following subsection, test the more general 

3- parameter model given in Dchncn & McLaugh- 
lin (2005). 

4- 1-2. Dehnen-McLaughlin's anisotropic 3-parameter 

model 

Dehnen & McLaughlin (2005, their equa- 
tion 46b) present a theoretically-motivated, 3- 
paramctcr model such that [a, (3, 7] = [2(2 — 
/?o)/9, (31 - 2/? )/9, (7 + 10/3 )/9], and the term 
/3o reflects the central (r = 0) anisotropy — a 
measure of the tangential to radial velocity dis- 
persion 2 . Setting 7' = (7 + 10/?o)/9, we have 
[a,/3, 7 ] = [(3 - 70/5,(18 - tO/5,7 7 ], and their 
density model can be written as 



p(r) 



2 6 Ps 



(r/r.)-r'[l + (r/r a )(3-T')/5]£ 



(17) 



As shown in Figure 10, for three of the six 
cluster-sized halos, this model has the greatest 
residual scatter of the four, 3-parametcr models 
tested here. For another two of the six cluster- 
sized halos it has the second greatest residual scat- 
ter. This model is also unable to match the curva- 
ture in the halos of the cold collapse models (Fig- 
ure 10). However, it does provide very good fits to 
the galaxy-sized halos, and actually has the small- 
est residual scatter for three of these halos (Ta- 
ble 2). 

The shallowest, inner, negative logarithmic 
slope of this model occurs when (3q = 0, giving 
a value of 7/9 « 0.78. For non-zero values of @o, 
this slope steepens roughly linearly with /3 - 

4.2. Sersic/Einasto model 

Sersic (1963, 1968) generalized de Vaucouleurs' 
(1948) i? 1 / 4 luminosity profile model by replac- 
ing the exponent 1/4 with 1/n, such that n was 
a free parameter that measured the 'shape' of a 
galaxy's luminosity profile. Using the observers' 
notion of 'concentration' (see the review in Gra- 
ham, Trujillo, & Caon 2001), the quantity n is 



2 Note: the quantities (3 and /3rj are not as related as their 
notation suggests. The former is the outermost, negative 
logarithmic slope of the density profile while the latter is 
the velocity anisotropy parameter at r = 0. 
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Fig. 9. — Residual profiles from application of the 
2-parameter (4/9, 31/9, 7/9) model (equation 16) 
from Dchncn & McLaughlin (2005, their equa- 
tion 20b) to our dark matter density profiles. 
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Fig. 10. — Residual profiles from application of 
the 3-parameter [(3 — 7') /5, (18 — 70/5, 7'] model 
(equation 17) from Dehnen & McLaughlin (2005, 
their equation 46b) to our dark matter density 
profiles. 
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monotonically related to how centrally concen- 
trated a galaxy's light profile is. With R denoting 
the projected radius, Sersic's R}l n model is often 
written as 

I(R) = I c exp {-&„ [(R/R e )^ n - l] } , (18) 

where I e is the (projected) intensity at the (pro- 
jected) effective radius R c . The term b n is not a 
parameter but a function of n and defined in such 
a way that R e encloses half of the (projected) to- 
tal galaxy light (Caon et al. 1993; see also Ciotti 
1991, his equation (1)). A good approximation 
when n > 0.5 is given in Prugniel & Simien (1997) 
as 

b n w In- 1/3 + 0.009876/n. (19) 

Assorted expressions related to the R}l n model 
can be found in Graham & Driver's (2005) review 
article. 

Despite the success of this model in describing 
the luminosity profiles of elliptical galaxies (e.g., 
Phillipps et al. 1998; Caon et al. 1993; D'Onofrio 
ct al. 1994; Young & Curric 1995; Graham et al. 
1996; Graham & Guzman 2003, and references 
therein) , it is nonetheless an empirical fitting func- 
tion with no commonly recognized physical basis. 
We are therefore free to explore the suitability of 
this function for describing the mass density pro- 
files, p{r) 1 of dark matter halos. Indeed, Einasto 
(1965, eq. 4; 1968, eq. 1.7; 1969, eq. 3.1) indepen- 
dently developed the functional form of Sersic's 
equation and used it to describe density profiles. 
More recent application of this profile to the mod- 
elling of density profiles can be found in Einasto 
& Haud (1989, their eq. 14) and Tenjes, Haud, 
& Einasto (1994, their eq. Al). Most recently, 
the same model has been applied by Navarro et 
al. (2004) and Merritt ct al. (2005) to character- 
ize dark matter halos, and Aceves, Velazquez, & 
Cruz's (2006) used it to describe merger remants 
in simulated disk galaxy collisions. 

To avoid potential confusion with Sersic's R}l n 
model, we define the following expression as 
"Einasto's r 1 /™ model" : 

p(r) = Pc exp {-d n [{r/r c f' n - l] } , (20) 

where r is the spatial (i.e., not projected) radius. 
The term d n , defined below, is a function of n such 
that p e is the density at the radius r c which defines 



a volume containing half of the total mass. The 
central density is finite and given by p(r = 0) = 
A=e d ". 

The integral of equation (20) over some volume 
gives the enclosed mass 3 , which is also finite and 
equal to 

M(r) =4ir p(f)f 2 df. (21) 
Jo 

This can be solved by using the substitution x = 
^(f/re) 1 /™ to give 



M(r) = Ai:nr 6 c p c e d "d n - M 1 (3n, x), (22) 

where j(3n,x) is the incomplete gamma function 
defined by 



7(3n,x) = f 
Jo 



(23) 



Replacing 7(3n, x) with T(3n) in equation (22) 
gives the total mass M tot . 

The value of d n , which we first saw in equa- 
tion (20), is obtained by solving T(3n) = 2 x 
7(3n, d n ), where F is the (complete) gamma func- 
tion. The value of d n can be well approximated 
(Mamon 2005, priv. comm.) by the expression 

d n w in - 1/3 + 0.0079/n, for n > 0.5 (24) 

(see Figure 11). 

In Paper II we recast Einasto's r 1 /" model using 
the radius r_ 2 , where the logarithmic slope of the 
density profile equals —2. 

Einasto's r 1 /™ model (see Einasto & Haud 
1989) was used in Navarro et al. (2004; their equa- 
tion 5) to fit their simulated dark matter halos. 
They obtained n w 1/(0.172 ± 0.032) « 6 ± 1.1. 
Subsequently, Merritt et al. (2005) showed that 
Einasto's r 1 /" model performed as well as the (1, 
3, 7) model, and gave better fits for the dwarf- 
and galaxy-sized halos, obtaining n w 5.6 ± 0.7. 
For a sample of galaxy-sized halos, Prada et al. 
(2005) obtained similar values of 6 — 7.5. 

Figure 12 shows the application of equation (20) 
to the A^-body halos of Section 3. A comparison 
with the (1, 3, 7) model fits in Figure 5 reveals 
that Einasto's model provides a better description 
for five of the six cluster-sized halos, three of the 
four galaxy-sized halos, and both of the spherical 
collapse halos. 



3 A similar expression is given in Mamon & Lokas 2005, their 
equation (A2); and Cardonc ct al. 2005, their equation (11). 
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for d n from equation (20), such that T(3n) = 
2j(3n,d n ), and the two approximations inset in 
the Figure. 
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Fig. 12. — Residual profiles from application of 
Einasto's r 1 /" model (equation 20) to our dark 
matter density profiles. 



Navarro et al. (2004) wrote "adjusting the pa- 
rameter [n] allows the profile to be tailored to each 
individual halo, resulting in improved fits" 4 . Such 
a breaking of structural homology (see Graham & 
Colless 1997 for an analogy with projected lumi- 
nosity profiles) replaces the notion that a universal 
density profile may exist. 

A number of useful expressions pertaining to 
Einasto's model, when used as a density profile 
(equation 20), are given in Cardone et al. (2005) 
and Mamon & Lokas (2005). In particular, Car- 
done et al. provide the gravitational potential, 
as well as approximations to the surface density 



and space velocity dispersion of the Einasto 



„l/n 

model, while Mamon & Lokas give approximations 
for the concentration parameter, central density, 
and M v i r i a i/M to tai- The nature of the inner pro- 
file slope of Einasto's r 1 /™ model and several other 
useful quantities are presented in Paper II. 



4.3. Prugniel-Simien model: 
jected Sersic R 1 / 11 model 



A depro- 



Merritt et al. (2005) tested how well a depro- 
jected Sersic i? 1 /" model fit p(r) from the Navarro 
et al. (2004) A-body halos. This was essentially 
the same as comparing the halo surface densities 
with Sersic's i? 1 /™ law. Prugniel & Simien (1997) 
presented a simple, analytical approximation to 
the deprojected Sersic law (their eq. B6): 



tie 



exp 



with 



P 



M 

T 



Le b b n{1 - p) 



-b(r/R e ) 



l/r, 



2R e T(n(3 - p))' 



(25) 



(26) 



Equation (25) is a generalization of equation (2) 
in Mellier & Mathez (1987), who considered only 
approximations to the deprojected i? 1 / 4 law. Mel- 
lier & Mathcz's model was itself a modification of 
equation (33) from Young (1976), which derived 
from the work of Povcda, Iturriaga, & Orozco 
(1960). 

In these expressions, i? e , n and b are understood 
to be essentially the same quantities that appear 



4 The value of n, equal to 1/a in Navarro et al.'s (2004) 
notation, ranged from 4.6 to 8.2 (Navarro et al. 2004, their 
table 3). 
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in the Sersic R}l n law that describes the projected 
density (equation 18). In fact, since equation (25) 
is not exactly a deprojected Sersic profile, the cor- 
respondence between the parameters will not be 
perfect. We follow the practice of earlier authors 
and define b to have the same relation to n as in 
equation (19). (For clarity, we have dropped the 
subscript n from b n .) Although the parameter p' is 
obtained from fitting the density profile, it can be 
defined in such a way that the total (finite) mass 
from equation (25) equals that from equation (18), 
giving equation (26). (We stress that the n in the 
Prugnicl-Simicn profile is not equivalent to the n 
in equation (20), Einasto's model.) 

This leaves the parameter p. We define p, like 
b, uniquely in terms of n: 

p= 1.0 - 0.6097/n + 0.05463/n 2 . (27) 

Lima Neto et al. (1999) derived this expression by 
requiring the projection of equation (25) to ap- 
proximate as closely as possible to the Sersic pro- 
file with the same (R e ,n), for 0.6 < n < 10 and 
10~ 2 < R/R e < 10 3 . 5 The accuracy of Prug- 
niel & Simien's (1997) approximation, using equa- 
tion (27) for p(n), is shown in Figure 13. 

Terzic & Graham (2005) give simple expres- 
sions, in terms of elementary and special func- 
tions, for the gravitational potential and force of a 
galaxy obeying the Prugniel-Simien law, and de- 
rive the spatial and line-of-sight velocity disper- 
sion profiles. 

One could also allow p to be a free parame- 
ter, creating a density profile that has any de- 
sired inner slope. For instance, setting p = 0, 
the Prugniel-Simien model reduces to the Einasto 
model. We do not explore that idea further here. 

The density at r = R c is given by p c — p'c~ b , 
while the projected surface density at R = R e , de- 
noted by I e , can be solved for using equation (26). 
Thus, one can immediately construct (a good ap- 
proximation to) the projected mass distribution, 
which will have a Sersic form (equation 18) with 
parameters (R e ,I e , and n). This allows the halo 
parameters to be directly compared with those of 
Sersic fits to luminous galaxies, which we do in Pa- 
per III. In Paper II we recast this model using the 

5 The value of p given in equation (27) is preferable to the 
value 1.0 - 0.6097/n + 0.05563/n 2 given in Marquez et al. 
(2000), (Lima Neto 2005, priv.comm.). 



radius where the logarithmic slope of the density 
profile equals —2. 

The mass profile (Terzic & Graham 2005, their 
Appendix A; see also Lima Neto et al. 1999 and 
Marquez et al. 2001), can be written as 

M(r) = 47r (0 'i? c 3 n&™ (p - 3) 7 (n(3 - p),Z) , (28) 

where Z = fe(r/i? ) 1 /" and 7(01,2;) is the incom- 
plete gamma function given in equation (23). The 
total mass is obtained by replacing 7 (n(3 — p), Z) 
with r (n(3 — p)), and the circular velocity is given 
by ?w (r) = sjGM{r)/r. 

In Figure 14, equation (25) has been applied to 
our dark matter profiles. The average (± standard 
deviation) of the shape parameter for the galaxy- 
sized and cluster-sized halos is n — 3.59(±0.65) 
and n = 2.89(±0.49), respectively. Merritt et al. 
(2005, their Table 1) found values of 3.40±0.36 
and 2.99±0.49 for their sample of galaxies and 
clusters, respectively, in good agreement with the 
results obtained here using a different set of N- 
body simulations and equation (25), rather than 
a numerically deprojected i? 1 /" light profile. 

Figure 14 reveals that CDM halos resemble 
galaxies (Merritt et al. 2005), since the projec- 
tion of the Prugniel-Simien model closely matches 
the Sersic R}l n model, and the latter is a good 
approximation to the luminosity profiles of stellar 
spheroids. Subject to vertical and horizontal scal- 
ing, CDM halos have similar mass distributions to 
elliptical galaxies with an absolute _B-band mag- 
nitude around — 18 ± 1 mag; these galaxies have 
n ~ 3 (see Graham & Guzman 2003, their fig- 
ure 9). This result was obscured until recently 
due to the use of different empirical models by ob- 
servers and modelers. 

Before moving on, we again remark that we 
have not explored potential refinements to the ex- 
pression (27) for the quantity p, but note that this 
could result in a better matching of the model to 
the simulated profiles at small radii. As the reso- 
lution of A-body clustering simulations continues 
to improve, it will make sense to explore such gen- 
eralizations. 

5. Model comparison: Which did best? 

Table 4 summarizes how well each parametric 
model performed by listing, for each type of halo, 
the rms value of A (equation 12a) for each set of 
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halos, given by 
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Fig. 13. — Logarithmic difference between the 
exact deprojection of Sersic's i? 1 /" model (equa- 
tion 18) and the approximation given by Prugniel 
& Simien (1997) in equation (25), using the val- 
ues of p and 6 given in equations (27) and (19), 
respectively. 
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Fig. 14. — Residual profiles from application of 
the Prugniel & Simien model (equation 25) to our 
dark matter density profiles. 



1 N 

N Z— 1 



(29) 



with N — 6, 4 and 2 for the cluster-sized, galaxy- 
sized, and spherical-collapse halos, respectively. A 
detailed description of each model's performance 
follows. 

The bowl- and hump-shaped residual profiles 
associated with the 2-parameter model of Burk- 
ert (1995) reveal this model's inability to describe 
the radial mass distribution in our simulated dark 
matter halos. The 2-parameter model of Dchncn 
& McLaughlin (2005) performs considerably bet- 
ter, although it too fails to describe the cold col- 
lapse systems and two of the six cluster-sized ha- 
los, specifically C09 and F09. Although this (4/9, 
31/9, 7/9) model never provides the best fit, it 
does equal or out-perform the NFW-likc (1, 3, 7) 
model in describing 3 of the 12 halos (A09, D12, 
G03). 

In general, all of the 3-parameter models per- 
form well (0.015 < A < 0.04 dcx) at fitting the 
A-body (non-collapse) halos. However, neither 
the (1, 3, 7) model nor the 3-parameter Dehnen- 
McLaughlin model can match the curvature in the 
density profiles of the cold collapse systems (Mil 
& M35). On the other hand, both Einasto's r 1 /" 
model and that from Prugniel & Simien give rea- 
sonably good fits (A ~ 0.05 dex) for these two 
halos. 

The Prugnicl-Simicn model provided the best 
overall description of the cluster-sized, iV-body ha- 
los. The (1, 3, 7) model and the 3-parameter 
Dehnen-McLaughlin model provided the best fit 
for only one cluster-sized, iV-body halo each, and 
even then the (1, 3, 7) model only just out- 
performed the Prugnicl-Simien model which gave 
the best fit for four of the six cluster-sized ha- 
los. For two of these halos, the size of the residual 
about the optimal Prugnicl-Simien fit was roughly 
half of the value obtained when using the (1, 3, 7) 
model. 

The implication of this result is that Sersic's 
i? 1 /™ model will describe the projected surface 
density of the cluster-sized, dark matter halos. 
Intriguingly, Demarco et al. (2003) and Durret, 
Lima Neto & Forman (2005) have observed that 
the (projected) hot X-ray gas distribution in clus- 
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Fig. 15. — Mass versus profile shape (1/n). For 
the galaxies and galaxy clusters, the shape pa- 
rameters n have come from the best-fitting Scrsic 
R 1 / 11 model to the (projected) luminosity- and 
X-ray profiles, respectively. The galaxy stellar 
masses, and cluster gas masses are shown here. 
For DM halos, the virial masses are shwon and 
the shape parameters have come from the best- 
fitting Prugnicl-Simicn model. (Note: The value 
of 1/n from the Prugniel-Simien model applied to 
a density profile is equivalent to the value of n 
from Sersic's model applied to the projected dis- 
tribution.) We are plotting baryonic properties for 
the galaxies alongside dark matter properties for 
the simulated halos. Filled stars: TV-body, dark 
matter halos from this paper; open plus signs: 
galaxy clusters from Demarco et al. (2003); dots: 
dwarf Elliptical (dE) galaxies from Binggcli & Jer- 
jen (1998); triangles: dE galaxies from Stiavclli 
et al. (2001); open stars: dE galaxies from Gra- 
ham & Guzman (2003); asterisk: intermediate to 
bright elliptical galaxies from Caon et al. (1993) 
and D'Onofrio et al. (1994). 



ters can indeed be described with Sersic's i? 1 /™ 
model; although the gas can at times display a 
rather unrelaxed behavior (Statler & Diehl 2006). 
Studies of gravitational lcnsing may therefore ben- 
efit from the use of Sersic's i? 1 /™ model for which 
the lensing equation has been solved (Cardone et 
al. 2004) and for which numerous other properties 
have previously been computed (Graham & Driver 
2005). 

With regard to the galaxy-sized, TV-body ha- 
los, the situation is somewhat different. Dchnen 
& McLaughlin's (2005) anisotropic 3-parameter 
model provided the best fit for three of the four 
profiles, with the Einasto r 1 /" model providing 
the best fit for the forth profile. We also observe 
that Einasto's model provided better fits than the 
(1, 3, 7) model for three of the four TV-body ha- 
los. If this observation holds, namely, that the 
Prugniel-Simien model describes the density pro- 
files of the cluster-sized halos best, while Dehnen 
& McLaughlin's 3-parameter model provides the 
best description of the galaxy-sized halos, it would 
imply that these halos do not have the same struc- 
tural form. Of course, even if the same model did 
provide the best fit for both types of halo, any 
variation in the value of the profile shape n, or 
central isotropy parameter (3 , would point toward 
the existence of nonhomology. 

While halos of different mass may be system- 
atically better described by different density laws, 
it is important to emphasize that a single den- 
sity law provides a good fit to all of the TV- 
body halos considered here. As Table 4 shows, 
Einasto's r 1 '" law has the smallest, or second- 
smallest, value of A rms for galaxy-sized, cluster- 
sized, and spherical-collapse halos. None of the 
other parametric models which we considered per- 
forms as well "across the board." The next best 
performer overall is the Prugniel-Simien profile. 

6. Discussion 

Figure 15 shows our TV-body halos, together 
with real elliptical galaxies and clusters, in the 
profile shape vs. mass plane. The profile shape 
parameter plotted there is either n from the Sersic 
R}l n model fit to the light profile, or the cor- 
responding parameter from the Prugniel-Simien 
model fit to the dark-matter density. Dynami- 
cal masses from the Demarco et al. (2003) study 
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of galaxy clusters are shown. We have also in- 
cluded the elliptical galaxy compilation in Gra- 
ham & Guzman (2003), converting their _B-band 
luminosities into solar masses using a stellar mass- 
to-light ratio of 5.3 (Worthey 1994, for a 12 Gyr 
old SSP), and an absolute £?-band magnitude for 
the Sun of 5.47 B-mag (Cox 2000). This ap- 
proach ignores the contribution from dark matter 
in galaxies. However, given the uncertainties on 
how Mt t/L varies with L (e.g., Trujillo, Burkert, 
& Bell 2004, and references therein) we prefer not 
to apply this correction, and note that the galaxy 
masses in Figure 15 only reflect the stellar mass. 

Figure 15 suggests that the simulated galaxy- 
sized halos have a different shape parameter, i.e. 
a different mass distribution, than the simulated 
cluster-sized halos. The same conclusion was 
reached by Merritt et al. (2005) who studied a 
different sample of TV-body halos. The sample 
of dwarf- and galaxy-sized halos from that paper 
had a mean (± standard deviation) 6 profile shape 
n = 3.04(±0.34), while the cluster-sized halos had 
n = 2.38(±0.25). We observe this same system- 
atic difference in our iV-body halos. Taking the 
profile shape n from the Prugniel-Simien model 
fits to the density profile (equivalent to the value 
of n obtained by fitting Sersic 's i? 1 /" model to the 
projected distribution) we find n = 3.59(±0.65) 
for our cluster-sized halos and n = 2.89(±0.49) for 
our galaxy-sized halos. A Student t test, without 
assuming equal variance in the two distributions, 
reveals the above means are different at the 88% 
level. Applying the same test to the data of Mer- 
ritt et al. (2005; their Table 1, column 2), which 
is double the size of our sample and also contains 
dwarf galaxy-sized halos, we find that the means 
arc different at the 99.98% level. We conclude that 
there is a significant mass dependence in the den- 
sity profiles of simulated dark-matter halos. Den- 
sity profiles of more massive halos exhibit more 
curvature (smaller n) on a log-log plot. 

The fact that n varies systematically with halo 
mass raises the question of which density scale 
and radial scale to use when characterizing halo 
structure. In the presence of a "universal" den- 
sity profile, the ratio between R e and r_2 (the 
radius where the logarithmic slope of the density 

'Reminder: the uncertainty on the mean is not equal to the 
standard deviation. 



profile equals —2, see Paper II) is a constant fac- 
tor, but with varying values of n this is not the 
case. This remark also holds for the scale den- 
sity, which is used to measure the contrast with 
the background density of the universe and pro- 
vides the so-called "halo concentration." This in 
turn raises the question of what "concentration" 
should actually be used, and whether systematic 
biases exist if one uses p-i rather than, say, p c . 
To reiterate this point: the density ratio between 
r = r_2 and r = R c depends on the profile shape 
n, and thus, apparently, on the halo mass. 

In Figure 16 we show how the use of r~2 and 
R c produce slightly different results in the size- 
density diagram (e.g., Figure 8 of Navarro et al. 
2004). The relation between size (or equivalently 
mass) and central concentration (or density) varies 
depending on how one chooses to measure the sizes 
of the halos. 

To better explore how the homology (i.e., uni- 
versality) of CDM halos is broken, it would be 
beneficial to analyze a large, low-resolution sam- 
ple of halos from a cosmological cube simulation 
in order to obtain good statistics. Moreover, the 
collective impact from differing degrees of virial- 
ization in the outer regions, possible debris wakes 
from larger structures, global ringing induced by 
the last major merger, triaxiality, and the presence 
of large subhalos could be quantified. 

7. Summary 

We presented a nonparametric algorithm for ex- 
tracting smooth and continuous representations of 
spherical density profiles from iV-body data, and 
applied it to a sample of simulated, dark matter 
halos. All halos exhibit a continuous variation of 
logarithmic density slope with radius; in the case 
of the ACDM halos, the variation of slope with 
radius is close to a power law. We then com- 
pared the ability of a variety of parametric models 
to reproduce the nonparametric /o(r)'s. Over the 
fitted radial range 0.01 < r/r v - lr < 1, both the 
Einasto r 1 /" model (identical in functional form 
to Sersic's model but expressed in terms of space, 
rather than projected, radius and density) and the 
Prugniel-Simien model (an analytical approxima- 
tion to a de-projected Sersic law) provide a better 
description of the data than the (1, 3, 7) model, 
i.e. the NFW-like double power-law model with in- 
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Table 4 

Residual scatter: rms values of A. 



Model 


Cluster-sized Galaxy-sized 
halos halos 


Spherical-collapse 
halos 




3-parameter models 






Einasto 


0.028 


0.026 


0.052 


Prugnicl-Simicn 


0.025 


0.030 


0.056 


(1,3,7) 


0.032 


0.028 


0.236 


Dchnen-McLaughlin (Eq.17) 


0.034 


0.023 


0.193 




2-parameter models 






Dehnen-McLaughlin (Eq.16) 


0.053 


0.032 


0.257 


NFW 


0.046 


0.046 


0.246 


Burkcrt 


0.218 


0.164 


0.217 



Note. — Col.(l): Model. Col. (2): rms of the 6 residual scatters, Arms (equation 29), 
for the cluster-sized halos. Col. (3): Similar to Col. (2) but for the 4 galaxy-sized halos. 
Col. (3): Similar to Col. (2) but for the 2 spherical collapse halos. For each halo type, 
the two models which perform the best arc highlighted in bold. 
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Fig. 16. — The density, p_2, where the logarith- 
mic slope of the density profile equals —2 is plot- 
ted against i) the radius where this occurs (open 
symbols), and ii) the effective radius (filled sym- 
bols) derived from the best-fitting Prugniel-Simien 
model (equation 25). Both p_ 2 and r_ 2 are also 
computed from the best-fitting Prugniel-Simien 
model, see Paper II. If a universal profile existed 
for these halos, then the vertical difference should 
be constant for all halos. 



ner slope 7. Moreover, unlike the (1, 3, 7) model, 
both of these models have finite total mass, and 
are also capable of describing the density profiles 
of halos formed from the cold collapse of a spher- 
ical over-density (Figure 6). 

The single function that provides the best over- 
all fit to the halo density profiles is Einasto's law, 
equation (20): 

p(r)=p c cxp[-d n [(r/r c ) 1 / n -l\} 

with d n defined as in equation (24). This con- 
clusion is consistent with that of an earlier study 
(Merritt et al. 2005) that was based on a different 
set of ./V-body halos. Typical values of the "shape" 
parameter n in equation (20) are 4 < n < 7 (Table 
1). Corresponding n values from Sersic profile fits 
to the projected (surface) density range from ~ 3 
to - 3.5 (Fig. 15). 

We propose that Einasto's model, equation (20), 
be more widely used to characterize the den- 
sity profiles of TV-body halos. As noted above, 
Einasto's model has already found application in 
a number of observationally-motivated studies of 
the distribution of mass in galaxies and galaxy 
clusters. We propose also that the suitability of 
Einasto's model for describing the luminous den- 
sity profiles of galaxies should be evaluated - ei- 
ther by projecting equation (20) onto the plane 
of the sky, or by comparing equation (20) directly 
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with deprojected luminosity profiles. Such a study 
could strengthen the already strong connection be- 
tween the density profiles of galaxies and TV-body 
dark-matter halos (Merritt et al. 2005). 

While equation (20) is a good description of 
all of the halo models considered here, we found 
that systematic differences do exist in the best- 
fit models that describe TV-body halos formed via 
hierarchical merging on the one hand, and those 
formed via spherical collapse on the other hand, in 
the sense that the latter have substantially smaller 
shape parameters, n w 3.3 (Table 1). That is, the 
density profiles in the cold collapse halos decline 
more quickly than r~ 3 at large radii, and have 
shallower inner profile slopes than those produced 
in simulations of hierarchical merging. 

With regard just to the non-collapse models, 
we also found systematic differences between the 
cluster- and galaxy-sized halos. The latter are 
slightly better fit by the 3-parameter Dehnen- 
McLaughlin model, and the former are slightly 
better fit by the Prugniel-Simien model (Table 
4). This, together with the observation that more 
massive halos tend to have smaller shape parame- 
ters n (Figure 15), suggests that there may not be 
a truly "universal" density profile that describes 
ACDM halos. 
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